!
!     PROGRAMME DIFEPH
!
!     Differences entre les ephemerides de la libration de la Lune :
!     Eph.DE245 - Sol.LLIB04
!     Eph.DE403 - Sol.LLIB04
!     Eph.DE405 - Sol.LLIB04
!
! --- Declarations -----------------------------------------------------
!
      implicit real*8 (a-h,o-z)
!
      character*50  fname
      dimension     v0(3),v1(3),v2(3),v3(3),dv(3,3)
!
! --- Parametres du calcul ---------------------------------------------
!
      dj0=2360400.5d0
      nbj=109000
      pas=20.d0
!
! --- Fichiers ---------------------------------------------------------
!
      nul1=11
      fname='LLIB04.DAT'
      open (nul1,file=fname,status='old')
      nul2=12
      fname='C:\TCHEBY\DE245LIB.AD'
      open (nul2,file=fname,status='old',access='direct',recl=968)
      nul3=13
      fname='C:\TCHEBY\DE403LIB.AD'
      open (nul3,file=fname,status='old',access='direct',recl=968)
      nul4=14
      fname='C:\TCHEBY\DE405LIB.AD'
      open (nul4,file=fname,status='old',access='direct',recl=968)
      nulout=20
      fname='DIFEPH.TXT'
      open (nulout,file=fname)
!
! --- Calcul des differences -------------------------------------------
!
      npas=nbj/pas
      dj=dj0
!
      do n=1,npas
         if (mod(n,10).eq.0) 
     &   write (*,'(2x,i5.5," / ",i5.5,2x,f9.1)') n,npas,dj
         call LLIB04 (dj,nul1,v0)
         call EPH245 (dj,nul2,v1)
         call EPH403 (dj,nul3,v2)
         call EPH405 (dj,nul4,v3)
         do i=1,3
            dv(i,1)=v1(i)-v0(i)
            dv(i,2)=v2(i)-v0(i)
            dv(i,3)=v3(i)-v0(i)
         enddo
         dv(3,2)=dv(3,2)+2.2d0
         dv(3,3)=dv(3,3)+3.9d0
         an=(dj-2451545.d0)/365.25d0+2000
         write (unit=nulout,'(f10.3,9f12.5)') an,dv
         dj=dj+pas
      enddo
!
      write (*,*) ' FINI '
!
      stop
      end
!
!
!
      subroutine LLIB04 (tjd,nul,var)
!-----------------------------------------------------------------------
!
!     Ref : JCGF 0407
!
! --- Object -----------------------------------------------------------
!
!     Computation of the lunar libration angles (solution LLIB04).
!
! --- Input ------------------------------------------------------------
!
!     tjd      Julian Date TDB (real*8).
!     nul      Logical unit number of the file LLIB04.DAT (integer).
!              This file has to be opened before the first call.
!
! --- Output -----------------------------------------------------------
!
!     var(3)   Libration angles : p1, p2, tau, in arcsecond (real*8).
!              Selenodetic system: axes of principal moments of inertia.              
!
! --- Declarations -----------------------------------------------------
!
      implicit real*8 (a-h,o-z)
!
      real*8,   parameter            :: cpi=3.141592653589793d0
      real*8,   parameter            :: rad=648000.d0/cpi
      real*8,   parameter            :: deg=cpi/180.d0
!
      logical                        :: xinit=.true.
!
      real*8,   dimension(3)         :: var
!
      real*8,   dimension(5)         :: t
      real*8,   dimension(2,500,3)   :: coef
      integer,  dimension(17,500,3)  :: iarg
      integer,  dimension(17)        :: iw
      integer,  dimension(3)         :: nterm
!
      real*8,   dimension(3,5)       :: w
      real*8,   dimension(5)         :: eart, peri
      real*8,   dimension(8,2)       :: p
      real*8,   dimension(4,5)       :: del
      real*8,   dimension(2)         :: zeta
      real*8,   dimension(3)         :: fli
!
      save
!
! --- Function : conversion between (deg,min,sec) and radians ----------
!
      DMS(ideg,imin,sec)=(ideg+imin/60.d0+sec/3600.d0)*deg
!
! --- Time TDT (unit: Julian century) ----------------------------------
!
      t(1) = 1.d0
      t(2) = (tjd-2451545.0d0)/36525.d0
      t(3) = t(2)*t(2)
      t(4) = t(3)*t(2)
      t(5) = t(4)*t(2)
!
! --- First call : initilization of lunar and earth-moon parameters ----
!
      if (xinit) then
!
         w(1,1)  = DMS(218,18,59.8782d0)  ! Mean mean longitude
         w(1,2)  = 1732559343.3328d0/rad  ! of the moon
         w(1,3)  =         -6.8700d0/rad
         w(1,4)  =   0.6604d-2/rad
         w(1,5)  =  -0.3169d-4/rad
!
         w(2,1)  = DMS( 83,21,11.6518d0)  ! Mean longitude 
         w(2,2)  =   14643420.3304d0/rad  ! of the lunar perigee
         w(2,3)  =        -38.2639d0/rad
         w(2,4)  =  -0.45047d-1/rad
         w(2,5)  =   0.21301d-3/rad
!
         w(3,1)  = DMS(125, 2,40.3265d0)  ! Mean longitude 
         w(3,2)  =   -6967919.8851d0/rad  ! of the lunar ascending node
         w(3,3)  =          6.3593d0/rad
         w(3,4)  =   0.7625d-2/rad
         w(3,5)  =  -0.3586d-4/rad
!
         eart(1) = DMS(100,27,59.1880d0)  ! Mean heliocentric
         eart(2) =  129597742.3016d0/rad  ! mean longitude 
         eart(3) =         -0.0202d0/rad  ! of the earth-moon barcenter
         eart(4) =   0.9d-5/rad
         eart(5) =   0.15d-6/rad
!
         peri(1) = DMS(102,56,14.4136d0)  ! Mean longitude
         peri(2) =       1161.2283d0/rad  ! of the perihelion
         peri(3) =          0.5327d0/rad  ! of the earth-moon barcenter
         peri(4) =  -0.138d-3/rad
         peri(5) =   0.d0
!
         fli(1)  =       44820417.d0/rad  ! Free libration
         fli(2)  =        1736493.d0/rad  ! arguments frequencies
         fli(3)  =    -5364715.227d0/rad
!
         p(1,1)  =   DMS(252,15, 3.25986d0)  ! Planetary longitudes J2000
         p(2,1)  =   DMS(181,58,47.28305d0)
         p(3,1)  =   eart(1)
         p(4,1)  =   DMS(355,25,59.78866d0)
         p(5,1)  =   DMS( 34,21, 5.34212d0)
         p(6,1)  =   DMS( 50, 4,38.89694d0)
         p(7,1)  =   DMS(314, 3,18.01841d0)
         p(8,1)  =   DMS(304,20,55.19575d0)
!
         p(1,2)  =   538101628.68898d0/rad  ! Planetary mean motions
         p(2,2)  =   210664136.43355d0/rad
         p(3,2)  =    eart(2)
         p(4,2)  =    68905077.59284d0/rad
         p(5,2)  =    10925660.42861d0/rad
         p(6,2)  =     4399609.65932d0/rad
         p(7,2)  =     1542481.19393d0/rad
         p(8,2)  =      786550.32074d0/rad
!
         dprec   =  0.d0
         preces  = (5029.0966d0+dprec)/rad  ! Constant of the precession
!
         do i=1,5
            del(4,i) = w(1,i)  - eart(i)    ! Arguments of Delaunay
            del(3,i) = w(1,i)  - w(3,i)
            del(1,i) = w(1,i)  - w(2,i)
            del(2,i) = eart(i) - peri(i)
         enddo
!
         del(4,1) = del(4,1)+cpi
         zeta(1)  = w(1,1)
         zeta(2)  = w(1,2)+preces
!
      endif
!
! --- Fist call : reading the lunar libration series LLIB04 ------------
!
      if (xinit) then
!
         read (nul,'(1x)')
         do iv=1,3
            read (nul,'(10x,i10)') nt
            nterm(iv)=nt
            do n=1,nt
               read (nul,'(5x,2d20.13,16i3,i5)') xs,xc,iw
               coef(1,n,iv)=xs; coef(2,n,iv)=xc
               do i=1,17
                  iarg(i,n,iv)=iw(i)
               enddo
            enddo
         enddo
!
         close (nul)
         xinit=.false.
!
      endif
!
! --- Computation of the libration angles p1, p2 and tau ---------------
!
      do iv=1,3
         nt=nterm(iv)
         var(iv)=0.d0
         do n=1,nt
            xs=coef(1,n,iv) ; xc=coef(2,n,iv)
            do i=1,17
               iw(i)=iarg(i,n,iv)
            enddo
            y=0.d0
            itest=iw(17)
            if (itest == 1 .or. itest == 3) then
               do i=1,4
                  do k=1,5
                     y=y+iw(i)*del(i,k)*t(k)
                  enddo
               enddo
            else
               do i=1,4
                  y=y+iw(i)*(del(i,1)+del(i,2)*t(2))
               enddo
            endif
            do i=5,12
               y=y+iw(i)*(p(i-4,1)+p(i-4,2)*t(2))
            enddo
            y=y+iw(13)*(zeta(1)+zeta(2)*t(2))
            do i=14,16
               y=y+iw(i)*fli(i-13)*t(2)
            enddo
            var(iv)=var(iv)+xs*sin(y)+xc*cos(y)
         enddo
      enddo
!
! --- Additive Poisson terms for tau (10-5-98) -------------------------
!
      y      = 18*p(2,1)-16*p(3,1)-del(1,1)+114.56550d0*deg
     &       + (18*p(2,2)-16*p(3,2)-del(1,2))*t(2)
      var(3) = var(3)+0.25d0*sin(y)*t(2)
      y      = del(2,1)+t(2)*del(2,2)
      var(3) = var(3)-0.23d0*sin(y)*t(2)
!
      return
      end subroutine LLIB04
!
!
!
      subroutine EPH245 (dj,nul,x)
!-----------------------------------------------------------------------
!
!     Ref : JCGF 0407
!
! --- Objet ------------------------------------------------------------
!
!     Calcul des angles de la libration de la lune (IN JPL).
!
! --- Entree -----------------------------------------------------------
!
!     dj       Date Julienne TDB (real*8).
!
! --- Sortie -----------------------------------------------------------
!
!     x(3)     Angles de la libration p1, p2, tau en " (real*8).
!              Repere: axes principaux d'inertie.
!
! --- Declarations -----------------------------------------------------
!
!
      implicit real*8 (a-h,o-z)
!
      real*8,  parameter       :: cpi=3.141592653589793d0
      real*8,  parameter       :: rad=648000/cpi
      real*8,  parameter       :: dxpi=2.d0*cpi
      real*8,  parameter       :: deg=cpi/180.d0
!
      real*8,  dimension(3)    :: x,y
      real*8,  dimension(5)    :: w
      real*8,  dimension(3,3)  :: rp,rq,a,b,rr
!
      DMS(ideg,imin,sec)=(ideg+imin/60.d0+sec/3600.d0)*deg
!
      w(1)  =  DMS(218,18,59.83482d0)
      w(2)  =  1732559343.35614d0/rad
      w(3)  =           -6.7996d0/rad
      w(4)  =           0.6604d-2/rad
      w(5)  =          -0.3169d-4/rad
!
      xphi =-0.07355d0
      xeps = 23*3600.d0+26*60.d0+21.40580d0
!
      gt = (dj-2451545.d0)/365250.d0
      t  = 10.d0*gt
!
!     Matrice de precession ecliptique entre J2000 et la date
!     pour le departure point (precession J. Laskar).
!
      p1 =  0.10180391d-4
      p2 =  0.47020439d-6
      p3 = -0.5417367d-9
      p4 = -0.2507948d-11
      p5 =  0.463486d-14
      q1 = -0.113469002d-3
      q2 =  0.12372674d-6
      q3 =  0.1265417d-8
      q4 = -0.1371808d-11
      q5 = -0.320334d-14
!
      pw   = ((((p5*t+p4)*t+p3)*t+p2)*t+p1)*t
      qw   = ((((q5*t+q4)*t+q3)*t+q2)*t+q1)*t
      pwqw = 2.d0*pw*qw
      pw2  = 1.d0-2.d0*pw*pw
      qw2  = 1.d0-2.d0*qw*qw
      ra   = 2.d0*SQRT(1.d0-pw*pw-qw*qw)
      pw   = pw*ra
      qw   = qw*ra
!
      rp(1,1) =  pw2
      rp(2,1) =  pwqw
      rp(3,1) =  pw
      rp(1,2) =  pwqw
      rp(2,2) =  qw2
      rp(3,2) = -qw
      rp(1,3) = -pw
      rp(2,3) =  qw
      rp(3,3) =  pw2+qw2-1
!
!     Matrice de passage du repere equatorial de l'IN
!     au repere dynamique ecliptique J2000
!
      phir    =  xphi/rad
      epsir   =  xeps/rad
      cphi    =  cos(phir)
      sphi    =  sin(phir)
      ceps    =  cos(epsir)
      seps    =  sin(epsir)
      rq(1,1) =  cphi
      rq(1,2) =  sphi
      rq(1,3) =  0.d0
      rq(2,1) = -ceps*sphi
      rq(2,2) =  ceps*cphi
      rq(2,3) =  seps
      rq(3,1) =  seps*sphi
      rq(3,2) = -seps*cphi
      rq(3,3) =  ceps
!
!     Calcul des angles d'Euler de la libration (IN JPL) 
!
      call TCHEB245 (dj,1,nul,y,ierr)
!
!     Matrice de passage des coordonnees equatoriales
!     du repere de l'IN aux coordonnees rapportees aux
!     axes principaux d'inertie de la lune d'apres l'IN
!
      phi    =  y(1)
      theta  =  y(2)
      psi    =  y(3)
      cphi   =  cos(phi)
      sphi   =  sin(phi)
      cpsi   =  cos(psi)
      spsi   =  sin(psi)
      cth    =  cos(theta)
      sth    =  sin(theta)
      a(1,1) =  cpsi*cphi-spsi*sphi*cth
      a(1,2) =  cpsi*sphi+spsi*cphi*cth
      a(1,3) =  spsi*sth
      a(2,1) = -spsi*cphi-cpsi*sphi*cth
      a(2,2) = -spsi*sphi+cpsi*cphi*cth
      a(2,3) =  cpsi*sth
      a(3,1) =  sphi*sth
      a(3,2) = -cphi*sth
      a(3,3) =  cth
!
!     Matrice de passage des coordonnees ecliptiques de la
!     date (departure point) aux coordonnees rapportees aux
!     axes principaux d'inertie de la lune d'apres l'IN
!
      do i=1,3
         do j=1,3
            d=0.d0
            do k=1,3
               d=d+rq(k,i)*rp(j,k)
            enddo
            b(i,j)=d
         enddo
      enddo
!
      do i=1,3
         do j=1,3
            d=0.d0
            do k=1,3
               d=d+a(i,k)*b(k,j)
            enddo
            rr(i,j)=d
         enddo
      enddo
!
!     Calcul des angles de la libration p1, p2, tau issus de l'IN
!
      p1   =   rr(1,3)*rad
      p2   =   rr(2,3)*rad
      d    =   1.d0+rr(3,3)
      stet =  -rr(2,1)+rr(2,3)*rr(3,1)/d
      ctet =   rr(1,1)-rr(1,3)*rr(3,1)/d
      s    =   atan2(stet,ctet)
      xl   =   (((w(5)*t+w(4))*t+w(3))*t+w(2))*t+w(1)
      xl   =   mod(xl,dxpi)
      s    =   s-xl-cpi
      s    =   mod(s,dxpi)
      if (s.gt.cpi)  s=s-dxpi
      if (s.lt.-cpi) s=s+dxpi
      tau=s*rad
!
      x(1)=p1 ; x(2)=p2 ; x(3)=tau
!
      return
      end subroutine EPH245
!
!
!
      subroutine EPH403 (dj,nul,x)
!-----------------------------------------------------------------------
!
!     Ref : JCGF 0407
!
! --- Objet ------------------------------------------------------------
!
!     Calcul des angles de la libration de la lune (IN JPL).
!
! --- Entree -----------------------------------------------------------
!
!     dj       Date Julienne TDB (real*8).
!
! --- Sortie -----------------------------------------------------------
!
!     x(3)     Angles de la libration p1, p2, tau en " (real*8).
!              Repere: axes principaux d'inertie.
!
! --- Declarations -----------------------------------------------------
!
!
      implicit real*8 (a-h,o-z)
!
      real*8,  parameter       :: cpi=3.141592653589793d0
      real*8,  parameter       :: rad=648000/cpi
      real*8,  parameter       :: dxpi=2.d0*cpi
      real*8,  parameter       :: deg=cpi/180.d0
!
      real*8,  dimension(3)    :: x,y
      real*8,  dimension(5)    :: w
      real*8,  dimension(3,3)  :: rp,rq,a,b,rr
!
      DMS(ideg,imin,sec)=(ideg+imin/60.d0+sec/3600.d0)*deg
!
      w(1)  =  DMS(218,18,59.87484d0)
      w(2)  =  1732559343.35624d0/rad
      w(3)  =           -6.7772d0/rad
      w(4)  =           0.6604d-2/rad
      w(5)  =          -0.3169d-4/rad
!
      xphi =-0.05294d0
      xeps = 23*3600.d0+26*60.d0+21.40928d0
!
      gt = (dj-2451545.d0)/365250.d0
      t  = 10.d0*gt
!
!     Matrice de precession ecliptique entre J2000 et la date
!     pour le departure point (precession J. Laskar).
!
      p1 =  0.10180391d-4
      p2 =  0.47020439d-6
      p3 = -0.5417367d-9
      p4 = -0.2507948d-11
      p5 =  0.463486d-14
      q1 = -0.113469002d-3
      q2 =  0.12372674d-6
      q3 =  0.1265417d-8
      q4 = -0.1371808d-11
      q5 = -0.320334d-14
!
      pw   = ((((p5*t+p4)*t+p3)*t+p2)*t+p1)*t
      qw   = ((((q5*t+q4)*t+q3)*t+q2)*t+q1)*t
      pwqw = 2.d0*pw*qw
      pw2  = 1.d0-2.d0*pw*pw
      qw2  = 1.d0-2.d0*qw*qw
      ra   = 2.d0*SQRT(1.d0-pw*pw-qw*qw)
      pw   = pw*ra
      qw   = qw*ra
!
      rp(1,1) =  pw2
      rp(2,1) =  pwqw
      rp(3,1) =  pw
      rp(1,2) =  pwqw
      rp(2,2) =  qw2
      rp(3,2) = -qw
      rp(1,3) = -pw
      rp(2,3) =  qw
      rp(3,3) =  pw2+qw2-1
!
!     Matrice de passage du repere equatorial de l'IN
!     au repere dynamique ecliptique J2000
!
      phir    =  xphi/rad
      epsir   =  xeps/rad
      cphi    =  cos(phir)
      sphi    =  sin(phir)
      ceps    =  cos(epsir)
      seps    =  sin(epsir)
      rq(1,1) =  cphi
      rq(1,2) =  sphi
      rq(1,3) =  0.d0
      rq(2,1) = -ceps*sphi
      rq(2,2) =  ceps*cphi
      rq(2,3) =  seps
      rq(3,1) =  seps*sphi
      rq(3,2) = -seps*cphi
      rq(3,3) =  ceps
!
!     Calcul des angles d'Euler de la libration (IN JPL) 
!
      call TCHEB403 (dj,1,nul,y,ierr)
!
!     Matrice de passage des coordonnees equatoriales
!     du repere de l'IN aux coordonnees rapportees aux
!     axes principaux d'inertie de la lune d'apres l'IN
!
      phi    =  y(1)
      theta  =  y(2)
      psi    =  y(3)
      cphi   =  cos(phi)
      sphi   =  sin(phi)
      cpsi   =  cos(psi)
      spsi   =  sin(psi)
      cth    =  cos(theta)
      sth    =  sin(theta)
      a(1,1) =  cpsi*cphi-spsi*sphi*cth
      a(1,2) =  cpsi*sphi+spsi*cphi*cth
      a(1,3) =  spsi*sth
      a(2,1) = -spsi*cphi-cpsi*sphi*cth
      a(2,2) = -spsi*sphi+cpsi*cphi*cth
      a(2,3) =  cpsi*sth
      a(3,1) =  sphi*sth
      a(3,2) = -cphi*sth
      a(3,3) =  cth
!
!     Matrice de passage des coordonnees ecliptiques de la
!     date (departure point) aux coordonnees rapportees aux
!     axes principaux d'inertie de la lune d'apres l'IN
!
      do i=1,3
         do j=1,3
            d=0.d0
            do k=1,3
               d=d+rq(k,i)*rp(j,k)
            enddo
            b(i,j)=d
         enddo
      enddo
!
      do i=1,3
         do j=1,3
            d=0.d0
            do k=1,3
               d=d+a(i,k)*b(k,j)
            enddo
            rr(i,j)=d
         enddo
      enddo
!
!     Calcul des angles de la libration p1, p2, tau issus de l'IN
!
      p1   =   rr(1,3)*rad
      p2   =   rr(2,3)*rad
      d    =   1.d0+rr(3,3)
      stet =  -rr(2,1)+rr(2,3)*rr(3,1)/d
      ctet =   rr(1,1)-rr(1,3)*rr(3,1)/d
      s    =   atan2(stet,ctet)
      xl   =   (((w(5)*t+w(4))*t+w(3))*t+w(2))*t+w(1)
      xl   =   mod(xl,dxpi)
      s    =   s-xl-cpi
      s    =   mod(s,dxpi)
      if (s.gt.cpi)  s=s-dxpi
      if (s.lt.-cpi) s=s+dxpi
      tau=s*rad
!
      x(1)=p1 ; x(2)=p2 ; x(3)=tau
!
      return
      end subroutine EPH403
!
!
!
      subroutine EPH405 (dj,nul,x)
!-----------------------------------------------------------------------
!
!     Ref : JCGF 0407
!
! --- Objet ------------------------------------------------------------
!
!     Calcul des angles de la libration de la lune (IN JPL).
!
! --- Entree -----------------------------------------------------------
!
!     dj       Date Julienne TDB (real*8).
!
! --- Sortie -----------------------------------------------------------
!
!     x(3)     Angles de la libration p1, p2, tau en " (real*8).
!              Repere: axes principaux d'inertie.
!
! --- Declarations -----------------------------------------------------
!
!
      implicit real*8 (a-h,o-z)
!
      real*8,  parameter       :: cpi=3.141592653589793d0
      real*8,  parameter       :: rad=648000/cpi
      real*8,  parameter       :: dxpi=2.d0*cpi
      real*8,  parameter       :: deg=cpi/180.d0
!
      real*8,  dimension(3)    :: x,y
      real*8,  dimension(5)    :: w
      real*8,  dimension(3,3)  :: rp,rq,a,b,rr
!
      DMS(ideg,imin,sec)=(ideg+imin/60.d0+sec/3600.d0)*deg
!
      w(1)  =  DMS(218,18,59.87267d0)
      w(2)  =  1732559343.32953d0/rad
      w(3)  =           -6.8368d0/rad
      w(4)  =           0.6604d-2/rad
      w(5)  =          -0.3169d-4/rad
!
      xphi =-0.05028d0
      xeps = 23*3600.d0+26*60.d0+21.40960d0
!
      gt = (dj-2451545.d0)/365250.d0
      t  = 10.d0*gt
!
!     Matrice de precession ecliptique entre J2000 et la date
!     pour le departure point (precession J. Laskar).
!
      p1 =  0.10180391d-4
      p2 =  0.47020439d-6
      p3 = -0.5417367d-9
      p4 = -0.2507948d-11
      p5 =  0.463486d-14
      q1 = -0.113469002d-3
      q2 =  0.12372674d-6
      q3 =  0.1265417d-8
      q4 = -0.1371808d-11
      q5 = -0.320334d-14
!
      pw   = ((((p5*t+p4)*t+p3)*t+p2)*t+p1)*t
      qw   = ((((q5*t+q4)*t+q3)*t+q2)*t+q1)*t
      pwqw = 2.d0*pw*qw
      pw2  = 1.d0-2.d0*pw*pw
      qw2  = 1.d0-2.d0*qw*qw
      ra   = 2.d0*SQRT(1.d0-pw*pw-qw*qw)
      pw   = pw*ra
      qw   = qw*ra
!
      rp(1,1) =  pw2
      rp(2,1) =  pwqw
      rp(3,1) =  pw
      rp(1,2) =  pwqw
      rp(2,2) =  qw2
      rp(3,2) = -qw
      rp(1,3) = -pw
      rp(2,3) =  qw
      rp(3,3) =  pw2+qw2-1
!
!     Matrice de passage du repere equatorial de l'IN
!     au repere dynamique ecliptique J2000
!
      phir    =  xphi/rad
      epsir   =  xeps/rad
      cphi    =  cos(phir)
      sphi    =  sin(phir)
      ceps    =  cos(epsir)
      seps    =  sin(epsir)
      rq(1,1) =  cphi
      rq(1,2) =  sphi
      rq(1,3) =  0.d0
      rq(2,1) = -ceps*sphi
      rq(2,2) =  ceps*cphi
      rq(2,3) =  seps
      rq(3,1) =  seps*sphi
      rq(3,2) = -seps*cphi
      rq(3,3) =  ceps
!
!     Calcul des angles d'Euler de la libration (IN JPL) 
!
      call TCHEB405 (dj,1,nul,y,ierr)
!
!     Matrice de passage des coordonnees equatoriales
!     du repere de l'IN aux coordonnees rapportees aux
!     axes principaux d'inertie de la lune d'apres l'IN
!
      phi    =  y(1)
      theta  =  y(2)
      psi    =  y(3)
      cphi   =  cos(phi)
      sphi   =  sin(phi)
      cpsi   =  cos(psi)
      spsi   =  sin(psi)
      cth    =  cos(theta)
      sth    =  sin(theta)
      a(1,1) =  cpsi*cphi-spsi*sphi*cth
      a(1,2) =  cpsi*sphi+spsi*cphi*cth
      a(1,3) =  spsi*sth
      a(2,1) = -spsi*cphi-cpsi*sphi*cth
      a(2,2) = -spsi*sphi+cpsi*cphi*cth
      a(2,3) =  cpsi*sth
      a(3,1) =  sphi*sth
      a(3,2) = -cphi*sth
      a(3,3) =  cth
!
!     Matrice de passage des coordonnees ecliptiques de la
!     date (departure point) aux coordonnees rapportees aux
!     axes principaux d'inertie de la lune d'apres l'IN
!
      do i=1,3
         do j=1,3
            d=0.d0
            do k=1,3
               d=d+rq(k,i)*rp(j,k)
            enddo
            b(i,j)=d
         enddo
      enddo
!
      do i=1,3
         do j=1,3
            d=0.d0
            do k=1,3
               d=d+a(i,k)*b(k,j)
            enddo
            rr(i,j)=d
         enddo
      enddo
!
!     Calcul des angles de la libration p1, p2, tau issus de l'IN
!
      p1   =   rr(1,3)*rad
      p2   =   rr(2,3)*rad
      d    =   1.d0+rr(3,3)
      stet =  -rr(2,1)+rr(2,3)*rr(3,1)/d
      ctet =   rr(1,1)-rr(1,3)*rr(3,1)/d
      s    =   atan2(stet,ctet)
      xl   =   (((w(5)*t+w(4))*t+w(3))*t+w(2))*t+w(1)
      xl   =   mod(xl,dxpi)
      s    =   s-xl-cpi
      s    =   mod(s,dxpi)
      if (s.gt.cpi)  s=s-dxpi
      if (s.lt.-cpi) s=s+dxpi
      tau=s*rad
!
      x(1)=p1 ; x(2)=p2 ; x(3)=tau
!
      return
      end subroutine EPH405
!
!
!
      subroutine TCHEB245 (tdj,ipv,nul,r,ierr)
!-----------------------------------------------------------------------
!
!     Ref : JCGF 0407
!
! --- Objet ------------------------------------------------------------
!
!     Calcul de la Libration de la Lune (angles d'Euler) a partir de
!     d'Ephemerides Tchebychev du JPL DE245.
!
! --- Entree -----------------------------------------------------------
!
!     tdj      Date julienne TDB (real*8).
!
!     ipv      Indice des coordonnees (entier).
!              ipv=1 positons.
!              ipv=2 positions et vitesses.
!
!     nul      Numero d'unite logique du fichier JPL DE245 (entier).
!
! --- Sortie -----------------------------------------------------------
!
!     r(6)     Table des resultats (real*8).
!              r(1) = phi (rad).
!              r(2) = theta (rad).
!              r(3) = psi (rad).
!              r(4) = vitesse  phi (rad/j).
!              r(5) = vitesse  theta (rad/j).
!              r(6) = vitesse  psi (rad/j).
!
!     ierr     Indice erreur (entier).
!              ierr=0  pas d'erreur.
!              ierr=10 erreur de fichier JPL DE245.
!              ierr=11 erreur dans tdj.
!              ierr=12 erreur dans ipv.
!
! --- Remarque ---------------------------------------------------------
!
!     Le fichier Ephemerides Tchebychev du JPL DE245 doit etre defini et
!     ouvert en lecture dans le moniteur (unite logique nul).
!
! --- Declarations -----------------------------------------------------
!
      implicit double precision (a-h,o-z)
      dimension r(6),resul(6),bloc(1000),coef(1000)
      dimension pc(15),vc(15)
      equivalence (coef(1),bloc(2))
!
! --- Initialisations --------------------------------------------------
!
      data ideb/0/,id0/0/
      data resul/6*0.d0/,tdj0/0.d0/,ipv0/0/,x0/2.d0/,xx0/100.d0/
      data pc/1.d0,0.d0,13*0.d0/
      data vc/0.d0,1.d0,13*0.d0/
      data dpi/6.283185307179586d0/
!
! --- Lecture du descriptif et test des parametres ---------------------
!
      if (ideb.eq.0) then
         ideb=1
      read (nul,rec=1,err=100) idf,tdeb,tfin,delta,dt,nbloc,ndmot,ncf
      idf=idf
      nbloc=nbloc
         ierr=11
         if (tdj.lt.tdeb.or.tdj.ge.tfin) return
         ierr=12
         if (ipv.lt.1.or.ipv.gt.2) return
      endif
      ncd=3*ipv
      ni=delta/dt
!
! --- Lecture des coefficients de Tchebychev ---------------------------
!
      id=(tdj-tdeb)/delta+2
      if (tdj.eq.tfin) id=id-1
      if (id.ne.id0) then
         read (nul,rec=id,err=100) (bloc(i),i=1,ndmot)
         dj1=bloc(1)
         id0=id
      endif
!
! --- Debut boucle calcul ----------------------------------------------
!
      if (tdj.ne.tdj0.or.ipv.gt.ipv0) then
         tdj0=tdj
         ipv0=ipv
         tf=tdj-dj1
         int=tf/dt
         int=min(int,ni-1)
         l=3*ncf*int+1
!
! ------ Calcul des positions ------------------------------------------
!
         x=2.d0*(tf-int*dt)/dt-1.d0
         if (x.ne.x0) then
            x0=x
            pc(2)=x
            dx=x+x
            do i=3,ncf
               pc(i)=dx*pc(i-1)-pc(i-2)
            enddo
         endif
         do i=1,3
            resul(i)=0.d0
            do j=1,ncf
               jp=ncf-j+1
               jt=l+ncf*(i-1)+jp-1
               resul(i)=resul(i)+pc(jp)*coef(jt)
            enddo
         enddo
         resul(3)=dmod(resul(3),dpi)
         if (resul(3).lt.0.d0) resul(3)=resul(3)+dpi
!
! ------ Calcul des vitesses -------------------------------------------
!
         if (ipv.eq.2) then
            xx=dx+dx
            if (xx.ne.xx0) then
               xx0=xx
               vc(3)=xx
               do i=4,ncf
                  vc(i)=dx*vc(i-1)+2.d0*pc(i-1)-vc(i-2)
               enddo
            endif
            do i=1,3
               resul(i+3)=0.d0
               do j=2,ncf
                  jp=ncf-j+2
                  jt=l+ncf*(i-1)+jp-1
                  resul(i+3)=resul(i+3)+vc(jp)*coef(jt)
               enddo
               resul(i+3)=resul(i+3)*2.d0/dt
            enddo
         endif
!
! --- Fin boucle calcul ------------------------------------------------
!
      endif
      do i=1,ncd
         r(i)=resul(i)
      enddo
      ierr=0
      return
!
! --- Erreur fichier ---------------------------------------------------
!
100   continue
      ierr=10
      return
!
      end subroutine TCHEB245
!
!
!
      subroutine TCHEB403 (tdj,ipv,nul,r,ierr)
!-----------------------------------------------------------------------
!
!     Ref : JCGF 0407
!
! --- Objet ------------------------------------------------------------
!
!     Calcul de la Libration de la Lune (angles d'Euler) a partir de
!     d'Ephemerides Tchebychev du JPL DE403.
!
! --- Entree -----------------------------------------------------------
!
!     tdj      Date julienne TDB (real*8).
!
!     ipv      Indice des coordonnees (entier).
!              ipv=1 positons.
!              ipv=2 positions et vitesses.
!
!     nul      Numero d'unite logique du fichier JPL DE403 (entier).
!
! --- Sortie -----------------------------------------------------------
!
!     r(6)     Table des resultats (real*8).
!              r(1) = phi (rad).
!              r(2) = theta (rad).
!              r(3) = psi (rad).
!              r(4) = vitesse  phi (rad/j).
!              r(5) = vitesse  theta (rad/j).
!              r(6) = vitesse  psi (rad/j).
!
!     ierr     Indice erreur (entier).
!              ierr=0  pas d'erreur.
!              ierr=10 erreur de fichier JPL DE403.
!              ierr=11 erreur dans tdj.
!              ierr=12 erreur dans ipv.
!
! --- Remarque ---------------------------------------------------------
!
!     Le fichier Ephemerides Tchebychev du JPL DE403 doit etre defini et
!     ouvert en lecture dans le moniteur (unite logique nul).
!
! --- Declarations -----------------------------------------------------
!
      implicit double precision (a-h,o-z)
      dimension r(6),resul(6),bloc(1000),coef(1000)
      dimension pc(15),vc(15)
      equivalence (coef(1),bloc(2))
!
! --- Initialisations --------------------------------------------------
!
      data ideb/0/,id0/0/
      data resul/6*0.d0/,tdj0/0.d0/,ipv0/0/,x0/2.d0/,xx0/100.d0/
      data pc/1.d0,0.d0,13*0.d0/
      data vc/0.d0,1.d0,13*0.d0/
      data dpi/6.283185307179586d0/
!
! --- Lecture du descriptif et test des parametres ---------------------
!
      if (ideb.eq.0) then
         ideb=1
      read (nul,rec=1,err=100) idf,tdeb,tfin,delta,dt,nbloc,ndmot,ncf
      idf=idf
      nbloc=nbloc
         ierr=11
         if (tdj.lt.tdeb.or.tdj.ge.tfin) return
         ierr=12
         if (ipv.lt.1.or.ipv.gt.2) return
      endif
      ncd=3*ipv
      ni=delta/dt
!
! --- Lecture des coefficients de Tchebychev ---------------------------
!
      id=(tdj-tdeb)/delta+2
      if (tdj.eq.tfin) id=id-1
      if (id.ne.id0) then
         read (nul,rec=id,err=100) (bloc(i),i=1,ndmot)
         dj1=bloc(1)
         id0=id
      endif
!
! --- Debut boucle calcul ----------------------------------------------
!
      if (tdj.ne.tdj0.or.ipv.gt.ipv0) then
         tdj0=tdj
         ipv0=ipv
         tf=tdj-dj1
         int=tf/dt
         int=min(int,ni-1)
         l=3*ncf*int+1
!
! ------ Calcul des positions ------------------------------------------
!
         x=2.d0*(tf-int*dt)/dt-1.d0
         if (x.ne.x0) then
            x0=x
            pc(2)=x
            dx=x+x
            do i=3,ncf
               pc(i)=dx*pc(i-1)-pc(i-2)
            enddo
         endif
         do i=1,3
            resul(i)=0.d0
            do j=1,ncf
               jp=ncf-j+1
               jt=l+ncf*(i-1)+jp-1
               resul(i)=resul(i)+pc(jp)*coef(jt)
            enddo
         enddo
         resul(3)=dmod(resul(3),dpi)
         if (resul(3).lt.0.d0) resul(3)=resul(3)+dpi
!
! ------ Calcul des vitesses -------------------------------------------
!
         if (ipv.eq.2) then
            xx=dx+dx
            if (xx.ne.xx0) then
               xx0=xx
               vc(3)=xx
               do i=4,ncf
                  vc(i)=dx*vc(i-1)+2.d0*pc(i-1)-vc(i-2)
               enddo
            endif
            do i=1,3
               resul(i+3)=0.d0
               do j=2,ncf
                  jp=ncf-j+2
                  jt=l+ncf*(i-1)+jp-1
                  resul(i+3)=resul(i+3)+vc(jp)*coef(jt)
               enddo
               resul(i+3)=resul(i+3)*2.d0/dt
            enddo
         endif
!
! --- Fin boucle calcul ------------------------------------------------
!
      endif
      do i=1,ncd
         r(i)=resul(i)
      enddo
      ierr=0
      return
!
! --- Erreur fichier ---------------------------------------------------
!
100   continue
      ierr=10
      return
!
      end subroutine TCHEB403
!
!
!
      subroutine TCHEB405 (tdj,ipv,nul,r,ierr)
!-----------------------------------------------------------------------
!
!     Ref : JCGF 0407
!
! --- Objet ------------------------------------------------------------
!
!     Calcul de la Libration de la Lune (angles d'Euler) a partir de
!     d'Ephemerides Tchebychev du JPL DE405.
!
! --- Entree -----------------------------------------------------------
!
!     tdj      Date julienne TDB (real*8).
!
!     ipv      Indice des coordonnees (entier).
!              ipv=1 positons.
!              ipv=2 positions et vitesses.
!
!     nul      Numero d'unite logique du fichier JPL DE405 (entier).
!
! --- Sortie -----------------------------------------------------------
!
!     r(6)     Table des resultats (real*8).
!              r(1) = phi (rad).
!              r(2) = theta (rad).
!              r(3) = psi (rad).
!              r(4) = vitesse  phi (rad/j).
!              r(5) = vitesse  theta (rad/j).
!              r(6) = vitesse  psi (rad/j).
!
!     ierr     Indice erreur (entier).
!              ierr=0  pas d'erreur.
!              ierr=10 erreur de fichier JPL DE405.
!              ierr=11 erreur dans tdj.
!              ierr=12 erreur dans ipv.
!
! --- Remarque ---------------------------------------------------------
!
!     Le fichier Ephemerides Tchebychev du JPL DE405 doit etre defini et
!     ouvert en lecture dans le moniteur (unite logique nul).
!
! --- Declarations -----------------------------------------------------
!
      implicit double precision (a-h,o-z)
      dimension r(6),resul(6),bloc(1000),coef(1000)
      dimension pc(15),vc(15)
      equivalence (coef(1),bloc(2))
!
! --- Initialisations --------------------------------------------------
!
      data ideb/0/,id0/0/
      data resul/6*0.d0/,tdj0/0.d0/,ipv0/0/,x0/2.d0/,xx0/100.d0/
      data pc/1.d0,0.d0,13*0.d0/
      data vc/0.d0,1.d0,13*0.d0/
      data dpi/6.283185307179586d0/
!
! --- Lecture du descriptif et test des parametres ---------------------
!
      if (ideb.eq.0) then
         ideb=1
      read (nul,rec=1,err=100) idf,tdeb,tfin,delta,dt,nbloc,ndmot,ncf
      idf=idf
      nbloc=nbloc
         ierr=11
         if (tdj.lt.tdeb.or.tdj.ge.tfin) return
         ierr=12
         if (ipv.lt.1.or.ipv.gt.2) return
      endif
      ncd=3*ipv
      ni=delta/dt
!
! --- Lecture des coefficients de Tchebychev ---------------------------
!
      id=(tdj-tdeb)/delta+2
      if (tdj.eq.tfin) id=id-1
      if (id.ne.id0) then
         read (nul,rec=id,err=100) (bloc(i),i=1,ndmot)
         dj1=bloc(1)
         id0=id
      endif
!
! --- Debut boucle calcul ----------------------------------------------
!
      if (tdj.ne.tdj0.or.ipv.gt.ipv0) then
         tdj0=tdj
         ipv0=ipv
         tf=tdj-dj1
         int=tf/dt
         int=min(int,ni-1)
         l=3*ncf*int+1
!
! ------ Calcul des positions ------------------------------------------
!
         x=2.d0*(tf-int*dt)/dt-1.d0
         if (x.ne.x0) then
            x0=x
            pc(2)=x
            dx=x+x
            do i=3,ncf
               pc(i)=dx*pc(i-1)-pc(i-2)
            enddo
         endif
         do i=1,3
            resul(i)=0.d0
            do j=1,ncf
               jp=ncf-j+1
               jt=l+ncf*(i-1)+jp-1
               resul(i)=resul(i)+pc(jp)*coef(jt)
            enddo
         enddo
         resul(3)=dmod(resul(3),dpi)
         if (resul(3).lt.0.d0) resul(3)=resul(3)+dpi
!
! ------ Calcul des vitesses -------------------------------------------
!
         if (ipv.eq.2) then
            xx=dx+dx
            if (xx.ne.xx0) then
               xx0=xx
               vc(3)=xx
               do i=4,ncf
                  vc(i)=dx*vc(i-1)+2.d0*pc(i-1)-vc(i-2)
               enddo
            endif
            do i=1,3
               resul(i+3)=0.d0
               do j=2,ncf
                  jp=ncf-j+2
                  jt=l+ncf*(i-1)+jp-1
                  resul(i+3)=resul(i+3)+vc(jp)*coef(jt)
               enddo
               resul(i+3)=resul(i+3)*2.d0/dt
            enddo
         endif
!
! --- Fin boucle calcul ------------------------------------------------
!
      endif
      do i=1,ncd
         r(i)=resul(i)
      enddo
      ierr=0
      return
!
! --- Erreur fichier ---------------------------------------------------
!
100   continue
      ierr=10
      return
!
      end subroutine TCHEB405
